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SYNOPSIS 


A comparative  study  of  step-by-atep  methods  which  are  commonly  used  in 
the  numeric*!  analysis  of  the  dynamic  response  of  structures  is  presented.  The 
method  of  analysis  is  based  on  the  general  theory  of  the  calculus  of  difference 
equations  and  the  algebra  of  matrices . The  available  step-by-step  techniques 
discussed  are  cla-sified  into  three  groups; 

1.  Acceleration  methods, 

2.  Difference  equation  methods, 

V Numerical  solutions  of  differential  equations. 

Comparisons  have  been  made  between  the  available  techniques  with  r«Bpect  to  the 
accuracy  of  a single  step,  propagation  of  errot’s  after  a length  of  time,  limita- 
tions imposed  by  instability. and  lack  of  convergence,  time  consumption,  and  self- 
checking provisions  of  the  procedures . The  purpose  of  the  work  has  been  to 
determine  the  range  of  'Applicability  of  the  various  techniques. 


COMPARISON  OF  STEP-BY-STEP  METHODS  FOR  ANALYZING  THE 
DYNAMIC  RESPONSE  OF  STRUCTURES 


I . INTRODUCTION 

1.1  Summary 

This  dissertation  is  concerned  with  the  analysis  'of  step-by-atep  methods 
commonly  used  in  numerical  solutions  of  the  dynamic  response  of  structures . 
Rigorous  solutions  are  not  always  possible  for  structures  with  non-linear  charac- 
teristics under  dynamic  loads  such  as  wind,  impact,  blast,  earthquake  or  vibratory 
motions,  particularly  in  the  case  of  multi.- degree -of -freedom  systems  with  plastic 
resistance  or  with  a varying  elastic  behavior  as  a function  of  time . Consequently 
a numerical  approach  is  indispensable  for  such  conditions  and  step-by-step  studies 
of  motion  with  respect  to  time  are  extremely  useful. 

The  purpose  of  this  dissertation  is  to  study  the  accuracy  and  range  of 
applicability  of  various  step-by-step  techniques  now  available  and  frequently 
used  in  problems  of  dynamic  response  of  structures . These  step-by-step  procedures 
may  be  classified  for  convenience  of  discussion  into  three  groups.  1.  Accelera- 
tion methods;  2.  Difference  equation  methods;  and  3-  Numerical  soJ.ution  of 
differential  equations.  In  the  acceleration  methods  the  displacement  and  velocity 
at  the  end  of  a time  interval  are  each  expressed  in  terms  of  the  displacement  and 
velocity  at  the  beginning  of  the  time  interval,  together  with  the  accelerations 
which  occur  at  the  ends  of  the  interval,  a lav  of  variation  of  the  acceleration 
within  this  time  interval  being  assumed.  The  acceleration  is  in  turn  governed 
by  the  differential  equation  of  motion  and  the  problem  may  theD  be  solved  by  an 
algebraic:  solution  of  simultaneous  equations  or  by  cut-and-try  iterations. 

The  second  group  of  methods  involves  the  application  of  finite  differ- 
ence equations  which  are  obtained  from  the  given  differential  equations  of  motion. 


Displacement  at  each  successive  step  during  the  motion  can  he  derived  from  the 
displacements  previously  obtained  by  means  of  finite  difference  equations  - For 
multi-degree -of -freedom  structures  the  solution  may  be  accomplished  by  solving  . 
a set  of  linear  simultaneous  equations  or  by  inverting  a matrix. 

The  third  group  of  methods  includes  conventional  devices  developed  by 
- mathematicians  for  the  numerical  solution  of  various  differential  equations  and 
which  are  quite  general  in  application.  They  may  be  adopted  even  in  more  compli- 
cated problems  than  those. involved  in  the  equation  of  motion  which  we  usually 
encounter  . . 

The  analyse  the  characteristics  of  each  of  the  available  techniques' it 
is  best  to  obtain  beforehand  an  algebraic  equation ' representing  each  of  the  various 
techniques,  of  step-by-step  procedures,  than  to  compare  It  with  the  rigorous  solu- 
tion of  the  differential  equation  of  motion  and  investigate  its  propagation  of 
errors.  This  can  be  done  in  one  of  the  following  ways.  First,  it  is  possible  to  .... 
express  the  approximate  solution  in  the  form  of  a finite  difference  equation  in 
terms  of  displacements  and  find  its  complementary  and  particular  solutions  by  means 
of  the  calculus  Of  finite  difference.  If  the  approximate  procedure  is  readily 
given  in  a finite  difference  form,  no  work  is  necessary  in  transfoi’ming  the 
original  procedure  into  finite  difference  equations.  Secondly,  the  set  of  linear 
equations  used  in  the  approximate  technique  may  also  lie  expressed  in  a matrix 
form  such  that  a column  matrix  consisting  of  displacement,  velocity  and  accelera- 
tion at  the  end  of  the  time  interval  is  equal  to  the  product  of  a square  matrix 
into  a column  matrix  consisting  of  displacement,  velocity  and  acceleration  at  the 
beginning  of  the  time  interval.  When  the  procedure  is  successively  carried  out 
n times,  the  square  matrix  multiplies  itself  to  the  nth  power  and  shows  the  re- 
lation between  the  initial  and  final  conditions.  The  former  way  is  more  simple 
as  far  as  mathematics  is  concerned  but  reveals  directly  the  dynamic  response  only 
in  displacement,  while  the  latter,  though  involving  more  algebra,  gives  not  only 


the  displacement,  hut  velocity  and  acceleration  as  well  if  desired. 

The  dynamic  analysis  of  a structure  is  usually  based  on  the  following 
assumptions: 

1.  The  mass  of  the  structure  may  be  represented  by  a number  of  separate 
concentrated  masses  supported  by  a flexible  and  weightless  framework, 

' 2,  The  resistance-deflection  relationship  of  the  structure  can  be  deter- 

mined beforehand  over  the  whole  range  of  action,  and  the  time  history  of  displace- 
ment or  external  forces ‘ is  known. 

Without  loss  of  generalization,  the  present  analysis  has  been  confined  to 
a single-degree-of-fx-eedom  system.  Nevertheless  the  motion  of  more  complicated 
multi -degree -of -freedom  systems  can  be  considered  as  being  made  up  of  the  motion 
in  several  modes,  each  mode  acting  as' a single -degree -of -freedom.  ~ 

, Generally  speaking,  accuracy  may  be  attained  if  the  time  interval  .is 
sufficiently  small  while  too  large  en  interval  may  produce,  very  misleading  results. 
However,'  since  different  degrees,  of  accuracy  can  result  from  different  methods  of 
application,  the  choice  of  time  interval  depends  upon  the  accuracy  desired  and  the 
amount  of  work  required; 

Acceleration  methods  need  no  special  training  for  their  application  since 
they  are  based  on  fundamental  concepts,  but  these  methods  are  always  handicapped 
by  the  criteria  of  convergence  and  stability.  The  constant  acceleration  method^1)* 
is  objectionable  because  of  its  rapid  divergence  of  amplitude.  Timoshenko's 
modified  acceleration  method^)  gives  better  results  than  that  of  constant  accel- 
eration, yet  the  frequency  error  is  still  appreciable.  It  is,  however,  free  from 
stability  difficulties  and  has'  no  enlarging  or  diminishing  effect  of  the  velocity 
response.  Newmark's  linear  acceleration  method  (2)  has  better  agreement  in 
frequency,  but  overshoots  a little  in  amplitude  due  to  the  enlarged  velocity 
response.  Newmark's  parabolic  acceleration  method^2)  has  even  better  agreement 


* Numbers  in  parentheses  refer  to  items  in  the  Bibliography. 


in  frequency,  tut  unfortunately  its  amplitude  diverges  exponentially,  and  it  is 
therefore  of  less  value  for  a long  lapse  of  time  in  spite  of  its  accuracy  in  the 
first  cycle  of  vibration.  Newroark's  /3-methcd(3)  maybe  regarded  as  a genera- 
lized acceleration  method,  introducing  a row  parameter  fi  in  the  displacement 
equation  so  as  to  control  the  effect  of  acceleration.  With  £<=1/4,  it  is  identical 
with  Timoshenko's  modified  method.  With  ^<=1/6,  it  is  the  same  as  the  linear 
acceleration  method.  It  corresponds  to  the  difference  equation  method  adopted  by 
I,evy(5)  yhen  0,  and  to  that  given  by  Salvador i(6)  whan  1/12.  The  great 
advantage lof  this  generalization  is  that  it  permits  a~ convenient  choice  of  .the 
time,  interval  determined  by  the  convergence  -criterion  during  the  operation.  ' 

Iilfferenceequatioh  methods  also  have  Criteria  for  stability.  These  ' 
procedure#: are  not?. self. -checking . A little  more  time  economy  nay  be  gained  since 
only  the  displacement  is  necessary  for  the  computation  and.,  the  velocity  may  b8 
disregarded  in  each  step  thus  saving  time  in  calculations . As  stated  before,  the 
difference  equations  adopted  by  Levy  and  Salvador i may  be  considered  Sb  identical 
to  Nevmark/s  /S  -method  when  =»  0 and  1/12  respectively, ",  except  that  the 
initial  conditions  are  treated  differently.  Houboit’s  method^)  is  said  to  be  an. 
improvement  over  Levy's  method,  since  it  employs  a cubic  curve  of  displacement  for 
the  difference  equation,  yet  it  suffers  from  the  converging  characteristic  of  the 
amplitude  and  from  a large  error  in  period.  The  computed  amplitude  of  Sh  undamped . 
system  aa  computed  by  this  method  will  decay  rapidly  after  a few  cycles  of  vibra- 
tion even  when  a small  time  interval  is  used. 

The  accuracy  of  the  numerical  solutions  of  differential  equations 
developed  by  Ruler,  Ruuge  and  Kutta  in  discussed  in  many  books  and  papers. 

The  application  of  these  methods  to  linear  vibration  problems  is  somewhat  time- 
consuming  in  comparison  with  the  methods  above  mentioned  particularly  in  multi- 
degree-cf-freedom  systems.  Runge -Kutta' s method  has  an  advantage  for  general 
applicability  in  that  it  is  always  stable  and  it  has  the  proper  criterion  for 


critical  damping  in  viscous  damping  conditions. 

Comparisons  of  true  amplitude  and  reriod  with  'pseudo1  or  computed  ampli- 
tude and  period  in  each  method  of  numerical  solution  are  made  to  investigate  the 
effects  of  length  of  time  interval,  natural  frequency  of  the  structure  and  other 
parameters.  Additional  discussion  of  these  factors  is  presented  in  later  chapters. 


II.-  GENERAL  METHODS  OF  ANALYSIS 


6. 


2 .1  Calculus  of  Finite  Difference  Equations 

* 

Analysis  may  be  made  for  each  method  by  expressing  the  given  differential  ■ 
equation  of  motion,  combined  with  the  procedure  of  operation,  into  a difference 
equation  Then  the  properties  of  this  difference  equation  represent  the  character- 
istics of  the  corresponding  numerical  method.  In  the  second  group  of  .available 
■ techniques  .described- in  the  last  chapter,  finite  difference  equations  are  readily 
formed  from  the  differential -equation  of  motion  by  replacing  the  higher  orders  of  ' 
derivatives  by  central  difference  patterns.  In  the  first  and  third  groups  of 
.available  techniques  more  algebraic  work  is  required  to  convert'  the  equations  of 

motion  into  a difference  equation.  However,  the  equations  of  operation  prescrib- 

• . • . ■ *.  «.•'  . " " ■ -■  • • .• 

. ing  the  given  motion  can  always  be  expressed ’in  terms  of  displacements  and  veloci-  : 

ties  In  a linear  relation,  and  can  easily  be  put  in  a difference  equation  form. 

In  acceleration  methods  the  equations  of  operation  may  at  first  contain  some 
'acceleration  terms  but  one  can  soon  eliminate  them  since  the  final  acceleration  - 
itself  can  be  expressed  in  terms  of  displacement,  velocity  and  initial  accelera- 
tion. Thus  if  the  equation  of  motion  is  given  in  the  form 

y + 2rpy  + py  * F(t)  (2.1.1) 

where  p is  the  natural  frequency  of  vibration  and  r the  coefficient  of  viscous 
damping  in  terms  of  p,  it  is  possible  to  represent  the  numerical  procedure  by  a 
. / finite  difference  equation  in  the  form 

<*ty»H  + asyn.,  *= bF(tiHI)+lkF(t„)+  bjffc.,)  (2.1,2) 

. or,  in  the  case  of  the  parabolic  acceJeration  method  or  Houbolt's  method, 

f a»y»  f y*-i  + o+yn* 


(2.1-3) 


The  solutions  of  these  difference  equations  are 


yn  » o, x,n  v- 

(2.1.4) 

and 

y„  * <2.  + 

+ eJ*j” 

(2.1.5) 

respectively,  where  x^, 

xn,  and 

xj  are  the  rootB  of  the  equation 

A.  Xa  + A*  X 

+ Am  S 

0 

/ n *]  C\ 

f — 

W 

\ t-  t u.  • v#  / 

or 

A,  X3  f X1 

f e?j  x 

+ 6.4  = 0 

(2.1.7) 

corresponding  to  Eqs.  (2.1.2)  and  (2.1.3),  and  c^,  cg,  and  cj  are  constants 
determined  from  the  initial  conditions. 

If; the  roots  X]_  and  xg  are  conjugate  complex  roots,  the, response  of  the 
numerical  procedure  is  periodic  although  there  may  exist  errors  in  both  amplitude 
and  frequency.'  One  the  other  hand  when  .all  the  x roots  are  real,  the  solution, 
hecomes  aperiodic  and  unstable.  By  'stable'  we  mean  that  the  response  of  the 
numerical' solution  remains  periodic  and  without  fluctuation  or  rapid  divergence  in 
amplitude.  As  far  as  time  period  is  concerned,  the  observation  of  these  roots 
serves  therefore  as  a criterion  of  stability.  Divergence  of  amplitude  may  also  be 
regarded  as. kind  of  instability  and  it  will  be  shown  in  later  chapters  that  it  is 
due  to  the  presence  of  a factor  with  an  exponential  power  of  time  which  occurs  in 
the  general  equation  of  response.  If  the  factor  equals  one,  the  amplitude  w-m 
neither  diverge  nor  converge  and  1b  therefore  stable.  When  the  factor  is  larger 
than  one,  the  amplitude  diverges  with  a rate  which  depends  on  the  magnitude  of  the 
factor.  Slow  divergence  is  generally  acceptable  in  some  problems,  it  is  determined 
by  the  allowable  error  in  amplitude  and  not  by  the  criterion  of  stability. 

Particular  solutions  of  these  difference  equationc  may  be  obtained  by  the 
calculus  of  finite  differences  though  sometimes  this  ui&y  involve  difficulties  in 
more  complicated  forcing  functions  , However,  an  approximation  can  generally  be 
made  by  expressing  the  forcing  function  in  a power  series  or  a Fourier  expansion 
which  is  always  solvable  in  this  kiud  of  finite  difference  equations. 

The  finite  difference  equations  (2.1.2)  and  (2.1.3)  consist  entirely  of 
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displacement  terms  and  therefore  the  general  solution  shows  only  the  response  in 
displacement  of  the  structure  at  the  end  of  the  time  interval  due  to  the  displace- 
ment and  velocity  at  the  beginning  of  the  time  interval  and  also  due  to  the  exteribr 
forces  if  there  are  any.  In  order  to  bring  out  the  response  in  velocity  of  the 
structure,  another  set  of  difference  equation  containing  all  velocity  terms  must 
be  formed  from  the  fundamental  equations  of  the  numerical  solution,  such  as 

t = b,  f(t/w)  + b^F  (t^jr  bjF(tr)./)  (2.1.8) 

°r  ♦ »,y,  * h,F(b,)+  bjF(t,.,)+  tkHtr.,) 

Similarly,  the  finite  difference  equations  may  also  contain  only  accelera- 
tion terms  in'  the  form  of 

t h Gi/n-t  s b,F(tn+i)+ bfF(tn)  t F (tn-i)  (2.1.10) 

or  6,%i  + *1%  * %y».i  a t>, F(tf,H)+  bf F(tn)+  bj F(U.,)+  bkF(t»mt)  (2.i.n) 


if , the  response  in  acceleration  is  required. 

All  the  results  of  numerical  solutions  are  henceforth  to  be  compared  with 
the  exact  solution.  In  the  present  analysis  the  motion  of  a structure  which  is 
assumed  to  be  of  elastic,  behavior  is  prescribed  by  the  well-known  differential 
equation  (2.1.1)  whose  general  solution  is  known  to  be 

y = trpt(A  COSt/t’r*  pt  + 8 s/'n <*//-  r*  pt)+  yp  (2.1.12) 


where  yp  ia  the  particular  solution,  and  A and  B are  constants  determined  from  the 
initial  conditions. 

For  free  vibration,  F(t)  ==  0,  yp  = 0,  ana  the  equation  of  motion  becomes 

y + 2rpy  + pFy  r 0 . (2.1.13) 


and 


The  solution  is 

y * e"^( yQ  yt  + 

y “ s c*s  pt  - 


(7>  A .IS') 
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For  free  vibration  without  damping,  the  equation  of  motion  can  further  be 
simplified  to 

y + p*y  9 ° (3.1.16) 

and  Its  solution  is 

a 

y * y*  KX pt  + -—-Stnpt  (2.1.17) 


and  J * /#  **spt  - py»  sin pt 


(2.1.18) 


3 . 2 Algebra  of  Matrices 

This  is  applicable  to  the  first  and  third  groups  of  methods  provided  that 
the  displacement,  velocity  and  acceleration  at  the  end  of  any  time  interval  can  be 
expressed. in  a linear  form  in  terms  of  the  displacement,  velocity,  acceleration 
and  some  other  parameters  at  the  beginning  of  the  time  interval.  For  example, 
when  yr  « y6  + <£,/gy0  f ^>i  yt  ^ 


y,  * Any*  * <H */®  * + a» 


(2.2.1) 


y>  * *n  % * *ti  y,  + y6  * oh 

a matrix  form  representing  these  linear  simultaneous  equations  can  be  written  as 


y/i 

y< 

K 

y> 

l / j 


&u  d/g  BLjf  On 

On  du  on  Ah 

&J/  &il  d-il  Of4 

o o o I 


Jo 

y* 

* « 

y» 

L / J 


(2.2.2) 


, or,  in  more  abbreviated 


fora, 


[/<]”  MM. 


(2.2.5) 

When  the  numerical  solution  is  carried  through  n successive  steps  of  equal  time 
duration,  with  the  final  displacement  function  of  a previous  time  interval  be- 
coming the  initial  condition  of  a new  time  Interval,  the  matrix  [Aj  operates  on 
itself  n tiffins,  so  that 


M - M'  M • 


(2.2.4) 


i 
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The  matrix  [a]q  can  he  expanded  by  means  of  the  Cayley-Hamilton  theorem  and 
Sylvester's  theorem  as  soon  as  the  characteristic  roots  are  obtained. 

The  characteristic  roots  of  the  matrix  [a]|  not  only  gives  the  expansion 

/if  Fa"!*1,  hut  alfln  dpterm-i  tips  t.Vip  nf  a+.nV,^  1 1 + v no  In  + Vie  fMn4+o 

L J ' -----  -•  * - ~~ 

difference  method  described  in  the  preceding  article.  The  presenae  of  a pair  of 
conjugate  complex  characteristic  roots  signifies  stability  and  periodic  response 
of  the  numerical  solution  while  all  real  roots  indicate  that  the  response  is  a- 
periodic. 

The  method  may  become  very  tedious  in  the  case  of  forced  vibrations  since 
the  presence  of  more  than  two  characteristic  roots  in  the  matrix  will  add  too  much 
algebraic  work  t-o  the  simplification  process  The  method  of  analysis  by  finite 
difference . equations  io  preferrable  in  this  case. 
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III.  ANALYSIS  OF  AVAILABLE  TECHNIQUES 

3 . .1  Acceleration  Methods 

3.1,1  Constant  Acceleration  Method 

The  basic  assumption  of  this  method  is  that  the  acceleration  of  the  mass 
in  motion  remains-  unchanged  throughout  a small  time  interval  and  is  equal  in 
direction  and  magnitude  to  the  acceleration  at  the  starting  point  of  the  concern- 
ing time  interval.  The  assumption  is  a rough  one,  and  provides  a rapid  but  in- 
accurate procedure.  The  error  in  this  method  is  so  large  that  it  is  seldom  uaed* 
A slight  modification  and  a little  more  work  improve  the  results  considerably 
The  advantage  of  speed  of  operation  cannot  compensate  for  the  loss  in  accuracy. 

Let  us  consider  first  a single  mass  in  free  vibration  without  damping. 
Then  from  elementary  mechanics  ve  obtain  the  following  expressions: 

/«/  = y* f (3  l.i.i) 

y*»  = y*  + ( 

* y*  + ynh  + ynJjjL  (3. i.i. a) 

where  h denotes  the  time  interval. 

Now  the  differentia],  equation  of  motion  for  a body  in  free  vibration 
without  damping  is 

/ * ply  ■ 0 (2.1.16) 

from  which  the  relation 

yn  « -/*V*  (3.1.1  3) 

is  obtained.  Substituting  in  Eqs . (3. 1.1.1)  and  (3. 1.1. 2),  we  get 

yn+i  - " f^hyn  f y,i  (j.i.i.4) 

and  y<H,  = ('-  ^)yn  * hy„  . (5. 1.1.5) 

From  these  relations  of  displacement  und  velocity,  one  gets  a finite  difference 


12. 


equation  in  terms  of  displacements  corresponding  to  the  computed  results  from  the 


constant  acceleration  method. 


y*ti  - 2 (t-  *7r)yn  + (ft  s 0 


(3  1.1.6) 


The  solution  of  this  finite  difference  equation,  together  vith  Eq. 
(3. 1.1.3),  yields  the  general  equation  for  the  displacement  predicted  by  the 


approximate  solution; 


0*  %')*(» 


k.Mh* 


where  ^ m astS:h  ptj 


* arc  cos 


■■  Sin  rt/u) 

/- 


(3- l.l. 7) 


(3 .1.1.8) 


Comparing  this  vith  J5q.  (2.1.1?)  of  the  exact  solution,  it  is  obvious  that  when 
the  time  interval  h is  very  small,  this  approximate  method  approaches  the  exact 
solution  as  a limit.  Since  the. time  interval  is  not  zero,  there,  is  an  error  in 
the  procedure.  We  split  the  resulting  error  into  tvo  parts,  one  is  the  error  in 
frequency  or  period,  the  other  is  in  amplitude.  The  phase  angle  np  should  bo 
equal  to  pt  if  the  solution  is  exact.  In  other  words,  the  exact  value  of  m 
should  be  pt/n  or  ph.  Hence  ve  obtain  a relationship  between  the  pseudo  period 


of  the  numerical  solution  and  the  true  period  of  the  exact  solution,  i.e. 

“ ^5T  (3.1. 1.9) 

The  frequency  has  an  error  of  3 percent  for  ph  ■ 05  and  of  10  percent  for  ph  « 


The  error  in  amplitude  is  objectionably  large  since  the  computed  dis- 
placement is  eubjected  to  a magnifying  factor  (t+  ^ which  diverges 

rapidly  with  the  number  of  steps  of  operations  n.  This  source  of  error  1b  dominant 
although  there  also  sxist  some  other  errors  in  the  coefficients  of  yc  and  yQ. 

The  coefficient  of  yn  becomes  (ft  l~  a,  fan  ftM~)  instead  of  1 

' ' />-  -7T  /'  („&}* 

and  varies  as  a function  oi'  n and  ph.  The  coefficient  of  y0  becomes 

instead  of  l/p,  which  sIbo  shows  a rapid  divergence  of  amplitude-  Fig.  1 


illustrates  the  rapid  divergence  of  the  envelope  of  amplitude  for  a single  moving 
mass  subject  to  yQ=0  and  y0=p. 

The  criterion  of  stability  shovs  that  ph  should  be  less  than  V Any 
vaiue  01  ph  greater  than  4 yields  aperiodic  response  cf  displacement  and  velocity. 

No  criterion  of  convergence  exists  for  this  method  since  one  operation  is 
sufficient  for  each  step  since  no  repetition  or  trial  necessary, 

For  free  vibration  with  viscous  damping,  the  analysis  is  more  complicated 
since  it  involves  one  more  parameter  r,  the  damping  factor  of  the  motion.  The 


difference  equation  now  becomes 


y,„  - * O-rphJO * #)»h  = a 


and  ite  solution  is 


OUili  '.AUU  IS  / • 

V-re  fEMt’Mu  ««,«,  •ft-  ,C- 

The  ratio  of  pseudo  period  to  true  period  becomes 

& - Bhi/EH 


(3.1.1.10) 


(3.1.1.11) 


(3.1.1. 12) 


Comparison  of  amplitudes  maybe  made  1‘romEqs.  (3  1.1.11)  and  (2,1.14)  Figs.  2 
and  3 shows  the  comparisons  of  period  and  amplitude  for  different  damping  factor 


Two  kinds  of  comparison  in  amplitude  are  made  for  all  techniques  des- 
cribed in  this  thesis.  The  first  one  deals  with  the  magnitude  of  the  envelope 
of  vibration  at  pt  » 1,  regardless  of  the  magnifying  effect  of  the  sine  term  in 
the  general  expression.  In  other  words,  this  is  a comparison  of  tha  exponential 
factor  which  multiplies  the  solution.  The  purpose  of  this  comparison  is  to 
demonstrate  the  rate  of  divergence  or  convergence  and  then  to  Judge  its  applica- 
bility. ThiB  ratio  varies  exponentially  with  pt,  and  therefore  the  amplitude 
ratio  at  any  instart  cf  time  may  be  found  by  its  exponential  relation  with  the 


14. 


ratio  at  pt  » 1 

Another  comparison  deals  with  the  peak  amplitudes  in  the  first  cycle  of 
vibration  due  to  an  initial  velocity  yQ  p.  The  first  peak  amplitude  occurs  at 

— ft  / > r 2!.^  v Vii  t 4*  ■*  +■  rro  v 0+0  fho  f m >.  <a  irolna  H'-ISI0  f'  ® ^ 

tione  due  to  the  error  ia  period  which  therefore  plays  an  important  role  in  the. 
pseudo  peak  amplitudes.  This  kind  of  comparison  may  give  a better  picture  both  of 
the  actual  oscillatory  motion  and  of  the  pseudo  motion  derived  from  the  approxi- 
mate techniques . 

The  criterion  for  stability  of  the  constant  acceleration  method  is 
expressed  by  the  following  equation. 

fh*  + 12  rph  - ~ 0 ' (3,1.1.13) 

This  shows  that  when . ph  > (4*/tr*4i  ■-  br)  , the  computed,  displacement  of 
motion  is  aperiodic  which  1 s not  true  for  r less  than  1.  For  the  critical  damp- 
ing condition,  i.e.  r *«  1,  ph  should  be  made  equal  to  or  greater  than  0-9282  in 
order  to  procure  an  aperiodic  response. 

The  constant  acceleration  method  is  too  crude  in  accuracy  and  therefore 
not  much  used  in  practice.  It  is  only  accurace  to  the  second  order  of  h and  errom 
may  arise  from  the  third  power  of  the  time  interval  since 

y>  - (*-  ^;)y*  * (<-rf>b)hye  (3.1.1.14) 

3*1-2  Timoshenko's  Modified  Acceleration  Method 

This  is  an  iracroved  method  obtained  from  the  last  one  by  modifying  the 
acceleration  of  the  moving  body.  The  acceleration  here  is  assumed  constant 
throughout  the  time  interval  and  equal  to  the  average  cf  ius  initial  and  final 
values.  The  frleaeutary  equations  of  momion  are  therefore 

y„„  * At  i(A  + Mt)h 


(3  1*2.1) 


y*n  a y»+t  i (y»  **■..  yin-Jh  ° y*  * y*  & + jfc yn h* t £ /n+i  fo* 


It  should,  be  noted  that  the  above  equations  nay  not  be  consistent.  When  corn-tired 
with  the  differential  equation  of  free  vibration  without  damping,  these  equations 


yield 


0*  <$)}»,  = -p'h*  * ('-  pT)r* 

(/+  ^)/„h  - 0-  *ir)y>  + by* 


The  difi’erence  equation  in  terms  of  displacement  now  becomes 


y*+f  *■  (z--u4n 


)y»  + y*-,  = o 


and  its  solution  is 


y„  = yo  tes  ryu  .+  ^ sin  /yu. 


where 


, Ph 

/*-*  s»7  /'  pvr 


Hriccs—^pr- 
!+  r!T 


Similarly  it  can  also  be  shown  that 

y„  = ye  ccsyu  - py«  sin  ny*  (3- 1.2.0; 

This  approximate  solution  has  no  error  in  amplitude]  neither  the  initial 

displacement  nor  the  initial  velocity  produces  errors  which  would  affect  the 

magnitude  of  displacement  or  velocity  thereafter  computed.  {See  Fig.  4)  The 

only  error  arising  in.  this  method  ia  due  to  the  difference  of  phase-angles  or 

the  discrepancy  in  period  or  frequency  which  can  be  expressed  by  the  equation 
Ts  _ bh  f „ 


T an.  sin  -jy- 
aud  is  plotted  in  Fig- 


(3J  -2.9) 


Another  advantage  of  this  method  is  that  no  criterion  for  stability  need 
be  imposed.  The  length  of  the  time  interval  can  be  chosen  corresponding  to  the 
accuracy  desired.  Unfortunately,  the  error  in  period  is  so  large  that  even  a 
time  interval  of  about  1/6  of  the  natural  period  will  result  an  8 percent  error 
in  frequency. 

If  viscous  damping  is  considered,  one  can  obtain  the  following  equations 


1 

r 
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cf  motion, 


Otfph-t  ^)y„t,  = -P'hy*  * ('-rpk-  4V/» 

( u-rph-+  = {t+rph~  ^jr  )y»  + ^y„ 

and  the  difference  equation  of  displacement 

H ■$*?)>'-»  - (*-$p)y<> r (‘'i * 0 


with  its  solution 


*-m -tfhr'jh. 


where  . /bh^JT  ^ /-  tit* 


(3.1.2.10) 
(3.1.2  11) 

(3-1.2.12) 
(3.1.2.13) 
(3 •1*2.14) 


The  criterion  of  at, ability  taay  be  expressed  as  /—  ^ 0 . This 

expression  is  independent  of  the  natural  frequency  and  time  interval,  i.e,  there 
is  no  value  of  ph  which  affects  the  stability  of  response.  The  stability  criterion 
is  identical  with  the  criterion  of  critical  damping. 

By  comparing  Eq.  (3.1-2.13)  with  the  exact  solution,  Eq.  (2.1-14),  one 
can  observe  that  both  period  and  amplitude  errors  exist  in  the  approximate  solu- 
tion. The  ratios  of  pseudo  to  actual  values  of  periods  and  amplitudes  are  as 
follows:  ■ 


ascun 


(6i)  = j 

rtj 
« *• 

• — ✓ 

l A Ut  pt  at  ( 

. / * rph  + fL»  J 

\A  Utptit  * ( 777pr+  7**r)  ‘ * 

Curves  are  plotted  for  these  equations  in  Pigs  5 a*id  6. 


(5-1  2.15) 
(3.1.2  16) 


3.I.3  Newmarg's  Linear  and  Parabolic  Acceleration  Method* 

These  methods  are  based  on  the  assumption  that  the  variation  of  the 
acceleration  of  the  mass  in  motion  is,  respectively,  linear  or  parabolic.  The 
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velocity  is  determined  from  a definite  integral  of  acceleration  over  the  range  of 
time,  and  the  displacement  is,  in  turn,  an  integral  of  the  velocity.  The  elemen- 
tary equations  of  motion  are  therefore  consistent. 

The  linear  acceleration  method  is  first  considered  The  equations  of 
motion  are 

//it/  * /n  * £y»  + y y»+i 

i * h * **  h*  •* 

and  y»H  a /n  thy„  + g y„  +•$/**/ 

»•  . 

In  case  of  free  vibration  without  damping,  y_  * - p9y 
equations  can  be  simplified  to  the  form 


(3. 1.3.1) 

(3. 1.3. 2) 

The 


~ * hy„ 

O'  *£)*»  - -/*  1 0-  , 

from  which  the.  difference  equation 

- (2-/C^r)y«  * y->  = 0 

is  derived.  The  solution  of  this  difference  equation  is 

/a  . 

y*m  y»  <***?  + pyTejP  **  v* , 

vhere  » * «tsi*  tfVllJjjL 


a asc  CCS 


/ 1 


w ■ 


(3. 1.3. 3) 

(5.1.3^) 
(3. 1.3.5) 

(3.1.3- 6) 

(3. 1-3.7)' 


Similarly,  = - /,/£>.//-  + yB  Cm  ty*.  (3. 1.3.8) 

Malting  comparisons  with  the  exact  solution  as  before,  the  error  in  the 
approximate  solution  consists  of  two  parts:  'the  error  in  period  and  the  error  in 
amplitude.  The  error  in  period  is  expressed  by  the  ratio 


7s ph 

<z 


The  amplitude  error  in  this  case  is  constant  for  a given  time  interval,  it  dees 
not  vary  with  the  lapse  of  time  and  ia  solely  a function  of  the  initial  velocity 
of  the  mass  If  the  mass  starts  from  rest  no  amplitude  error  will  occur  Hence 
the  amplitude  error  depends  on  the  proportion  of  initial  displacement  and  velocity 
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The  error  will  be  large  for  a motion  due  to  a small  initial  displacement  and  large 
initial  velocity;  and  will  be  small  for  a motion  due  to  a small  initial  velocity. 
Fig.  7 shows  the  response  of  a single  mass  due  to  y^  -•  0*  and  yQ  = p. 

The  criterion  of  stability  ia  found  to  be  ph<^\/l2.  Any  value  of  ph 
greater  than  a/I2  will  result  in  aperiodic  response  by  the  approximate  solution. 

In  the  case  of  free  vibration  with  damping,  the  difference  equatior  of 
the  procedure  becomes 


The  solution  is 


where 


..  ..  .... 


(3.10.11) 


,yh>  ~**us ,/(%  eff.  Tyh1  ■ (j.i.3.12) 


The  criterion  for  stability  is  now 

/-  ?*-  * 0 for  ?<i  (5.13  13) 

Characteristics  of , the  pseudo  period  and  pseudo  amplitude  are  shown  on 
Figs , 8 and  9* 

The  parabolic  acceleration  method  differs  from  the  above  method  in  assum- 
ing that  the  acceleration  has  a parabolic  variation;  thuB  yH  » t k,h  t kthl  • 
The  procedure  of  operation  is  similar  to  the  linear  acceleration  method  except 
that  one  more  initial  condition  is  required  to  start  with.  That  is,  one  needs 
two  previously  known  steps  to  carry  out  a new  step.  The  procedure  may  be  started 
m one  of  the  following  ways* 

]. . Use  linear  acceleration  for  the  starting  interval. 

2 Use  linear  acceleration  for  a starting  interval  which  is  only  half 
ao  long  as  the  regular  interval,  then  got  a special  parabolic  acceleration  in- 
terval half  as  long  as  the  regular  interval  and  proceed. 
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The  equations  of  motion  from  elementary  mechanics  are 

/»  * 1-  <>  s'y„.,  i sy.)  (3-3  !■*) 

/y?  * /yw  -f  A»-'^  + fa  (- fag  I1  /tfj'ysy  +3pn)  (J.I.;  15) 


%***  J. 


The  difference  equation  corresponding  to  the  above  equations  is 

(»  f)y»  - (*--&#% * (<+£&')»  - - 9 • 


->  4. 

^ ry*V'-“ 


(y.  .1.5.16) 


The  solution  of  this  equation  can  he  expressed  in  the  form  of  Eq.  (2.1.5)  which 
can  also  he  written  as 

y„  * Ap*+  R”(8  cosn/A  + Cs/n  */*.)>  (3.1.5.17) 

where  A,  B and  G are  constants  determined  from  initial  conditions,'  p,  R and  u are 
functions  of  ph,  The  value  of  p:is  small  and  terms  containing  it  are  of  less  im- 
portance than  other  terms  of  the  expression.  R is  dominant  since  it  is  greater 
than  1 and  of  exponential  power  n as  steps  of  time  proceed.  Therefore,  the  pro- 
cedure is  of  a divergent  nature  and  is  not  desirable  over  a long  lapse  of  time 
Furthermore,  the  criterion  of  stability  ph  <[3.557758  limits  the  applicability  of 
long  time  intervals.  Although  the  parabolic  acceleration  method  is  more  accurate 
than  the  linear  one  as  far  as  the  first  cycle  of  oscillation  is  concerned,  the 
errors  in  amplitude  enlarge  rapidly  from  the  second  cycle  onward.  (See  Fig  10.) 
Figs.  11  and  12  show  how  Ig/T  and  kjk  vary  with  ph,  neglecting  the  first  term  of 
Eq.  (3.1.5.17). 


An  attempt  has  been  made  to  Improve  the  accuracy  of  the  linear  accelera- 
tion method  and  to  lessen  the  work  of  computation  of  the  parabolic  acceleration 
method  by  applying  linear  and  parabolic  accelerations  alternately  in  successive 
steps,  i.e.,  using  linear  acceleration  in  the  first  step,  parabolic  in  the  second, 
linear  in  the  third,  parabolic  in  the  fourth,  and  so  on.  The  result,  as  one  may 
expect,  turns  out  to  be  intermediate  between  the  two  methods . The  difference 
equation  for  the  displacements  at  an  eveo  number  of  steps  in  the  case  of  free 
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vibration  with  no  damping  is  written  as 


y it- &£t2y  > / * -Asfh *+  7ft ft* V -0 

y**"0  “ 7z gt£)(f+ mp)  y**  * t,  pjk)Y/*^ii)t  ‘-®h“  0 (31. 5.18) 


o T) 


and  its  solution  is 


.1.4  - 

/J  n % r >r  '.c  - 3 l>  n / * 


^ “ * fa  * y^/6  - TZStyb V ttoMpP-Mph*  Sm^] 


vhere  /.  = antm  ft, , /3MzMES53mE?3.  $2 
r V 2204+  1344 ph'i  388  f&i4+  S6p*h(+  ; 


' V 2$04+/344p*h'i388r»4+Sbp 

= ««■*«,  *Ca*-*ir»<l+sW) 


3 fW 


and  ^ — 


V £i<?4  * /J  44  368  pW+St^+jffh9 

V^Ja4  7 /J44ptht+  388  P*h  % J-g  s'/V  * 


(8+  p’h'XL+P'h3) 


(3.1.3.19) 


(3.1.3.20) 


(3.1.3.21) 


The  values  of  T6/T  and  As/A  are  plotted  against  ph  in  Figs.  11  and  12. 

It  is  evident  that  the  propagation  of ‘error  is  divergent.  The  criterion  of  stabi- 
lity of  this  alternate  linear-parabolic  acceleration  method  for  free  vibration 
without  damping  is  found  to  be  ph  * 1.6171  which  is  quite  unfavorable  for  long 
time  intervals. 


3.1.4  Newmark's  ^-Method 


This  is  a generalization  of  first  degree  acceleration  methods  obtained  by 
introducing  a parameter  £ into  the  elementary  equations  of  kinematics.  Thus 

y**,  ~ A + £ A *■  j y**>  (3. 1.4.1) 

y**t  • y«  + y*h  + (£-8)y«k*  + e'ymih1  (3. 1.4.2) 

It  is  obvious  here  that  this  iB  equivalent  to  Timoshenko's  modified 
.-method  when  p a 1/4  and  to  Newrtark 1 s linear  acceleration  method  -hen  (3  ■ 1/6. 

In  the  case  of  free  vibration  without  damping,  the  difference  equation 


of  displacements  appears  as 

a*/  ~ (*-  i+jtpw)?” r “ 15  • 


(3  1-4.3) 
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The  general  solution  of  the  difference  equation  becomes 

y„  - y,  usy*  V f ^=4=^.  s,n  T 


vter  a * **6Sm  > t Jtnt/.i = Ore car  — 7777 


/-(i-tip'h' 


(3. 1.4.4) 


(VI. 4. 5) 


By  comparing  the  pseudo  period  with  the  true  period,  it  can  be  seen  that 
the  ratio,  Tfi/T  may  be  expressed  in  a series  form  of 

_2_„  BL 

r a 

* /+£fa~)/6t-jj£c(7£opt-/?o{S  + /7)p4^----  (3.1  a.6) 


Consequently,  for  p * l/12,  the  ratio  T3/l  Will  be  closest  to  unity  for  any 
arbitrary  value  of  ph.  This  means  that  0 « 1/12  will  give  the  best  result  in 


period. 


The  computed  amplitude  ia  neither  divergent  nor  convergent  as  time  pro- 


gresses, although  some  error  is  involved  in  the  response  to  an  initial  velocity. 
Eq.  (3. 1.4, 4)  shown  that  the. term  containing  y0  does  not. contain  any  error  in 
amplitude,  but  the  one  with  y0  is  amplified  by  the  factor  Vi- 
which  is  only  dependent  on  p and  ph  and  is  free  from  influence  of  the  proceeding 
time.  Fig.  13  illustrates  the  variation  of  the  velocity  amplification  factor 


/-  (t-fl)p'i 


with  ph  for  different  values  of 


The  criterion  for  stability  is 

P'h'  , s 
Y+pp'k*  V.  * 


(3- 1.4.7) 


When  p » 1/4,  this  condition  ia  always  fulfilled  for  any  value  of  ph  and  therefore 
no  stability  criterion  need  be  imposed.  With  p^l/4,  no  ph  can  possibly  satisfy 
the  criterion.  With  p^l/4,  the  larger  0,  the  longer  the  time  interval  which  can 
be  used.  For  0 n 0,  the  critical  ph  is  equal  to  2. 

The  criterion  for  convergence  ia 


Ph  < ^/T  (3. 3. 4.8) 

which  shows  that  a larger  & permits  a smaller  applicable  time  interval.  For  3 « 0 


i 
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the  procedure  immediately  converges  regardiesa  of  the  time  interval  used. 

For  problems  of  free  vibration  with  viscous  damping,  the  difference  equa- 
tion of  the  procedure  may  be  vrittcn  as 

. , l \ f ' /•  rt*h  ' 

v+iijpwj fat " * (/-  ~ (5. 1.4. 9) 


Its  solution  is 
/ 

where 


on  is  • 

/t'rft*  fifty  \ '//-r'- (i-fi)pW  V (J*1*^* 


10) 


/tt-  SS  tircsm 


« arcus 


phj~- “IW&tL 

* Y0+ fi/Wf-  ryt>‘ 

VC'*  wwf-r'p'k' 


(J.  1.1* /ll) 

Now  the  error  in  period  is  not  o.’ily  dependent  on  p and  ph,  but  also  on  r. 

Again,  if  the  ratio  Ts/T.  is  expressed  in  a series  form  as 
jg_  m Pbs/j-r*' 

7 „ . A r/«> 3 - / j -+CW+ '>» t ICib'M  f (3.1.4.12) 

£4  f / - ry  f 

it  will  be  found  that  ^ s (l  + +t9  ~ &r*)/IZ  (l - 2 tz)  .(3. 1.1*. 13) 

is  the  best  value  as  far  us  period  is  concerned . 

The  error  in  amplitude  is  dominated  by  the  exponential  factor  multiplying 
the  whole  expression  in  Sq.  (3. 1.,1*. 10)  particularly  after  a long  sequence  of  time. 
Neglecting  the  coefficients  which  combine  with  y0  and  y0  in  the  expression,  one 
may  compare  amplitudes  by  taking  the  ratio 

(3.1.4.14) 

6 ( l + rph  t fiph'i 


(At) 

( A Jrfptsi 


For  best  agreement  in  amplitude, 


Jl 


The  criterion  for  stability  is 

pV  < -AQ^- 


(3-1.4  15) 

for  r^l.  Except  when  [5  =■  1 /h>  the  numerical  procedure  does  not  present  an  agree- 


ment on  the  criterion  of  critical  damping,  greater  discrepancy  usually  occutb  for 
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larger  time  Intervals.  Critical  damping  may  occur  in- the  numerical  procedure  even 
for  low  values  of  r if  too  long  a time  interval  is  taken. 


The  criterion  for  convergence  is 


rph  + fifth*  < / 


(3.1.4.161 


There  is  one  more  restriction  of  this  method  in  the  viscous  damping  pro* 

ffih 

blem.  A degenerate  case  of  the  difference  equation  will  occur  if  —j? 

is  equal  to  +}  or  -1.  Under  this  condition  the  method  is  not  workable.  However, 

when  - — rfi.h  — s / 

i+nrt‘ 


rph  + pp'h1  = zfipW  , 


ve  see  that  the  convergence  criterion  is  violated.  On  the  other  hand,  when 

r£h.  t . 

the  spring  constant  is  negative. 

The  3-method  is  also  applicable  to  forced  vibrations  as  represented  by 
Eq.  (2.1.1)  with  good  accuracy.  The  error  due  to  the  presence  of  a forcing 
function  F(t)  may  be  seen  by  comparison  with  the  exact  particular  solution.  Now, 
let  yp  be  the  particular  solution  of  the  given  differential  equation  of  motion, 
Eq.  (2.1.1),  and  let  yp  be  the  corresponding  solution  obtained  by  the  numerical 
procedure.  Analysis  is  made  for  an  undamped  system  of  single  mass  subject  to  an 


applied  force  i’(t)  as  follows: 

Given  the  equation  of  motion 

/ + ply  = F(t) 

the  exact  solution  is 


ys  (y*-yt>a)z>spt  + yP 


(3.1.4,17) 


(3,1.4.18) 


where  y is  the  particular  solution, 
P 


The  corresponding  difference  equation  when  using  the  p-method  is  found 


tc  be 


y*"  * (*-  - 7. ffa-i)] 


(3-1.4.19) 
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and  its  solution  is 


y»  a A c.osn/4  + 8 n/t  + y? 


where  yn  is  the  particular  solution  of  the  difference  equation  and 

/-nrrr^n  /_  /u. a) 

* « ^ - .«»«,  • 7 ‘rffrr 


a - y.  - »i 


and  8 «■ 


wP 


Comparing  Eqe.  (3.1.4.18)  and  (3.1.4.20),  we  see  that  when  h-^0,  if  yp-*>  J , the 
numerical  solution  will  approach  the  exact  solution.  The  following  comparisons 
are  made  for  different  kind  of  forcing  functions; 

(a)  The  forcing  function  is  a polynomial  in  time,  i.e, 
f(t)  S CL*  + d,t  + *%t*  t dji1* (3. 1.4. 2] 

then 

+ [ar^>+  S+p'1'- 

-ksJgJL  - 

>U-£&'+  T-t-gn 1-  *+ry *«+....]? 


+ f*  3‘2k,(ti.  Sd’Sz  focis  7'6‘5-4'2'2  7+ It 

fUj  ^ ^ J 

+ ^ 4>ikt,Gi{  i>  5'4029  t* 

* [a.  S-4  k, as,  7-6  Sjhat  ldJL**±h£t +....]£*  (3-1.4. 23; 


+ - 
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vhere 


6"  / ■/  isOtp-Opw 

ki  = i+  i ( & fi-  >)Plhl+ 7£o  (-360/3*-  top  + OP4t/* 

k+  - ! T i &(sopl-4o$+ 

kt~  t+i  (w-  >W+  jk  (<440^40^ 


/,  8/4,400 


(l&l4,46of-6o48ocfrss,92op*-/,<02ofi  * l)p‘/f9 


From  the  above  ve  see  that  the  numerical  solution  is  exact  for  a third 
degree  polynomial  at  any  value  of  $,  and  for  a fifth  degree  polynomial  when,  0 * 


(b)  The  forcing  function  is  trigonometric,,  i.e. 

f(t)  » A sm  at  + B tos  bt 


then  yp  - jrral  * n«-TT  601  ^ 


and  yp  z 


pK  a*  p 

A 4/h  at 

^ . 2 2 in* ah 


' (t+cos  ah-Z/i sin* ah  )h  * y O+cos 


B Cos  bt 

2 sin  * bh 


(c)  Exponential  forcing  functions,  as 

f(l)  = Of 

- *L 

* “ f>'+(bSa.)‘ 


» “ T* 


'-2a*  + l 


fe-2)a*  + /] 


(3*1.4.24) 


(3.1,4.25) 


(3.1.4.26) 


(3.1.4.87) 

(3.1.4.28) 


(3.1.4.29) 


In  all  cases  above,  vhen  b -*>  0,  yp-**  y^/  &hfi  therefore  the  numerical  solu- 
tion approaches  the  exact  solution  as  a limit. 

The  0 -method  may  also  be  applied  to  other  form  of  motion  than  the  one  re- 
presented by  Eq<  (2.1.1).  Consider  the  notion  prescribed  by  the  linear  differen- 


tial equation 


y - y *»  t 


(3.1.4.30) 
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the  exact  solution  of  which  is  known  to  he 

p = ce*  - (t  + i)  (5.1.4.31) 

where  c is  a constant  determined  from  the  initial  conditions.  Applying  the  0-' 
ctct-hcd  pc  wiiis  p2*c^l2!!i,  cns  Titids  "tins  dld,?srisnc,5  ^cuspio—  of  dpaviiAr^rcpnt.  to 


t 


* (i-4fi)h%  » $[ti  (/■  4f)h]t„., 


and  its  general  solution  is 

y * c( ttht-JL-pr)  -(t+0 

The  heat  value  of  0 for  this  case  is  1/6. 


(3.1.4.32) 


(3.1.4.35) 


3.2  Methods  of  Finite  Differences 
3 .2.1  Levy's  Method 

A method  proposed  by  Levy  replaces  the  second  derivative  yn  in  the  equa- 
tion of  motion  by  finite  differences,  ~jjf  (y#+t  " */»  + )/»•<)  . For  a free 

vibration  with  no  damping  prescribed  by  the  jiven  differential  aquation 

y + p*y  * * (2.1.16) 

by  substituting  yi  - -h  (y™  - * y*-> ) , one  obtains 

yn+t  - (2 - f>'hl)y*  * >-/  * 0 (3. 2. 1.1) 

which  is  obviously  identical  with  Navmark ' s generalized  acceleration  method  for  0*0, 
except  that  the  treatment  of  initial  velocity  is  different.  The  general  solution 
of  Eq.  (3.2  1)  Is 

X,  * Acesnyu  + & &'n  (3-2.12) 

where  ^ _ CUTC^'n  pb+/f-  = arc  COS  ( f - ) 

(3. 2- 1.3) 


and  A and  S are  constants  determined  from  initial  conditions,  with  ph  ^ 2 as 
stability  criterion.  There  is  no  difficulty  in  finding  A which  is  always  equal 
to  y„,  but  trouble  arises  in  the  determination  of  B which  depends  on  the 
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interpretation  of  the  initial  velocity,  or  in  other  words,  on  the  vay  in  which  y, 
is  obtained  from  the  known  values  of  yQ  and  y0.  When  yQ  n 0,  one  may  assume  that 
V-L  3 y.u  then  solve  simultaneously  with  Eq.  (J.2.1.1)  to  get  y^.  Otherwise,  if 


If  one  begins  with  Nevmark's  3-mathod  of  3 » Q for  the  first  step  in  order 


to  obtain  y^,  then 

y * 0-  + H 

and  the  result  will  be  the  same  as  Hewmark’s  3-method  for  3 «*  C,  i.e., 

3 - fiyjl  TjJE 

If,  taking  the  formula  of  elementary  mechanics 

y,  * y»  + hy0 

we  obtain  & + ^*r /• 


(3. 2. 1.4) 


(3. 2. 1.5) 


(3-C.1.5) 


(3. 2. 1.7) 


The  result  is,  of  course,  less  accurate. 
On.  taking  /#  = x * 


for  the  first  half  time  interval 


and  getting  y^  from  th9  difference  equation  of  the  half  time  interval 

Z%')i  * #> 

^ (j.2.1.8) 

This  result  is  generally  better  than  that  of  Eq.  (3. 2. 1.7). 

On  assuming  y_^  » 0 and  solving  for  y^, 

3 - * , (3. 2.1.9) 

This  is  only  true  when  y3  .s  \ .. ■ , and  is  only  used  when 

J n /y  /- 

is  uncertain. 

In  viscous  damping  problems,  by  replacing  yn  by  £%  ( y^,  " y»>i)  and 
fTi  by  jfi  ( /ah  **  2y*  + y*.,)  in  the  differential  equation  of  motion, 

y + 2 rpy  + F>*y  * o , (2.1.13) 

one  may  obtain  the  difference  equation 

0+ffh)y»,  - (z-pl»l)y*  * ('->Ph)yi-i  = « (5-24.10) 


(3. 2. 1.9) 

, and  is  only  used  when 


-v.  — 
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which  is  idervtial  to  Eql  (3. 1.4. 9)  for  p ■>  0.  The  result  will  be  the  same  as  that 
of  NevTPaeh  ■ it  generalized,  acceleration  method  for  p * 0 with  exception  of - the  treat- 
ment of  initial  velocity  The  discussion  is  therefore  not  repeated  here. 

Fux-  comparison  uf  pseudo  and  true  periods  and  amplitudes  of  this  method, 
see  Figs  . 14  and  13 . 


3.2.2  Galvadori's  Method 


This  is. on  application  of  a;  procedure  due  to  Fox  vibration  problems.  The 
second  derivative  y is  replaced  by  the  first  two  terms  of  the  central  difference 


expansion, 


7 s it)* 


(3.2  2.1) 


vher&.  A1  and  A*  are  the  second  and  fourth  central,  differences  of  y 

Axyn  * y*tt  - iyn  + y*-r 

* y*+i  - 4y»tt  + by*  ~ ly*.,  + /*.# 

f 4ft*  ) 

Operating  then  with  ( 1+  Jz  i bn  each  term  of  the  equation 

y + p'y  o 

Dropping  the  sixth-difference  terms,  the  equation  may  be  simplified  to 

>'  - (t ' fnw)  /» f >'  * o 


(3. 2. 2. 2) 


which  is  obviously  the  sane  bb  Nevmark's  generalized  acceleration  method  for 


P®  1/12.  (See  Eq.  (3*1. 4. 3).)  The  solution  is  in  the  form 

y„  .a  a ef/yu  + BSthrya 


/a.  a arc  sin 


°hS^~W 


it  & 


- arc  ers  — 


>+  £ 

1 * 71 


(3. 2. 2. 3) 


(3*2.?. 4) 


The  constants  A and  B are  determined  from  the  initial  conditions . For  free  vibra- 
tion^ A is  always  equal  to  y0,  while  B depends  on  the  way  in  which  y^  is  obtained 
from  the  initial  velocity  y0 . 

If  the  procedure  is  started  with  Nevmark ' s (3-method  for  p =>  1/12  'until 
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is  obtained,  that  ie  if 

X -77 


* ty]  , 


(5. 2. 2,5) 


the  result  will  be  the  same  as  ♦•hat  given  by  flevmark's  p-aethod  for  p ■ 1/12. 

y* 


a — 
a - 


•0.2.2. b) 


If  we  take  the  formula  of  elementary  mechanics 

y,  = y»  + > (3.2.2.?) 

B a — ~y/'mpjfir* — (3. 2.2.8) 

which  involves  mole  error  than  the  previous  result.  Hovevav,  the  accuracy  of  the 


then 


velocity  response  can  be  improved  if  the  initial  velocity  is  properly  treated. 

In  the  case  of  forced  vibrations  with  a forcing  function  F(t),  the 
difference  equation  of  Salvadorl's  method  beooaas 

y*«  * y*t  (3,2,2,9) 

which  is  again  the  Base  as  Newark's  method  for  P ■ 1/12.  (Bee  Eq.  3.1. h.,19).) 
This  method  is  accurate  to  the  order  of  h^  provided  that  the  motion  starts  from 
rest.  If  the  motion  begins  with  a finite  velocity,  the  treatment  of  the  initial 
velocity  for  the  difference  equation  determines  the  accuracy.  The  discussion  of 
this  method  is  included  in  Newark's  method  In  previous  and  later  chapters,  and 
is  not  restated  here. 

Salvadori  treats  the  damped  motion  problem  in  the  saae  way  as  Levy  does 
by  transforming  the  damped  motion  equation  into  an  algebraic  equation  by  the 
substitution  of  central  averaged  differences  for  the  derivatives.  The  difference 
equation  thus  formed  is 

(i+rph)yM-  t (t-rfik)y,.t  * 0 (3.2.2.10) 

which  is  the  same  as  Eq.  (3. 2.1. 9).  There  is  a difference  in;  procedure  when 
applied  to  a multi-degree-of -freedom  structure  but  it  does  not  affect  the  nature 


0?  the  errors . 
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3.2.3  Eoubolt 1 s Method  ~ 

Houbolt's  method  is  based  on  the  assumption  of  a cubic  curve  for  the  dis- 
placement of  the  moving  body,  considering  that  four  successive  ordinates  can  be 
passed  through  by  a cubuc  curve.  With  this  assumption  the  following  diff ei-cuce 


equations  for  the  final  derivative  may  be  obtained 

y*  * "fa  ( 2yn  - sy»-i  * 4y»-*  - /*• a)  (3. 2. 3.1) 

> s th( //>  “ 16 >' + 3 y*'!  “ (3. 2. 3. 2) 

The  derivatives  at  the  third  of  the  four  ordinates  are  sometimes  of  use 
and  are  as  follows:- 

y»  - -'*/»  + y™)  ■ (3-2.30) 

y*  ~ + iy*  - ty»-i  + /*••)  0.2.3.10, 

For  free  vibration  without  damping,  substitute  Eq  (3.2. 3. 1)  into  the 
differential  equation  of  motion 

y + f>ly  m 0 (2.1.16) 

The  following  difference  equation  is  obtained. 

(t+f>xhl)yn  - Syn.,  + 4/r.t- y*-3  m a (3. 2.3.5)' 

1 ’■  , 

Its  solution  is 

> * - c,x!*  + c,xS  ftf/V  (3. 2. 3.6). 

where  x.p  x^}  and  are  the  roots  of  the  equation 

(z  + phl)x*  - Sx*  + 4X  - / * 0 (3. 2.3-7) 


It  can  be  shown  hy  the  theory  of  equations  that  this  equation  contains  one  real 
root  and  two  conjugate  complex  roots  for  any  value  of  ph.  Therefore  the  solution 
has  always  an  oscillating  nature  and  no  criterion  of  stability  governs  the  choice 
of  the  time  interval,  although  the  amplitude  may  be  damped  out  very  rapidly  ae 
time  proceeds  (see  Fig.  18).  Eq.  (3. 2. 3 .6)  may  also  be  written  in  the  form 

> « At'***  f t-*rr(8  cos  cpt  + C swept)  (3 ,2,3.8) 

where  a,  b,  and  c are  all  functions  of  ph.  Here  A,  E and  C are  constants  deter- 
mine.- from  the  initial  conditions  The  first  term  of  the  equation  is  negligibly 


small  while  the  second  term  multiplied  by  e * ' is  dominant,  Table  (3. 2=3.1)  lists 
the  numerical  values  of  a,  b and  c for  various  ph,  Fig.  19  shows  the  ratio  of 
pseudo  period  to  true  period  T /T,  and  Fig.  20  shows  the  ratio  of  amplitudes  Aa/A. 

The  disadvantages  of  this  method  are  two-fold.  First,  It  needs  one  more 
initial  condition  to  start  with.  Although  this  can  be  found  by  taking  account  of 
the  initial  acceleration  it  also  involves  effort  to  trace  the  hack  differences 
including  the  simultaneous  solution  of  equations  Eqs.  (3. 2. 3.1)  to  (3, 2,3.4).  If 
the  initial ' conditions  are  awkwardly  treated,  e.g.  by  making  tho  ‘assumption  that 
the  fictitious  displacements  at  t » »h  and  t s -2h  are  zero,  the  error  introduced 
by  these  erroneous  assumptions  would  be  greater  than  that  of  the  method  itself. 

The  step-by-step  evaluation  of  succeeding  displacements  cannot  proceed  In  a 
straightforward,  manner  until  three  initial  displacements  have  been  established. 
Secondly,  the  amplitude  of  a slightly  damped  system  decreases  so  rapidly  even  for 
a time  interval  of  about  one-sixth  of  the  natural  period,  that  the  amplitude  is 
reduced  50  percent  after  one  and  one-hajf  cycluB  of  vibration.  (See  Fig.  18.) 

Finer  time  intervals  and  thus  more  computational  effort  must  be  uBed  to  reduce 
the  damping  effect  of  the  procedure. 

For  free  vibration  with  vIbcous  damping,,  the  difference  equation  of  this 
method  becomes 

(2+j^rp(i+p*h)y„  - (s+  & <*ph)yn.,  + ( 4+ 3rph)y*.z-(h  jrph)y*.3  « 0 (3. 2. 3. 9) 

The  solution  is  in  the  form  of  Eq.  (3. 2.3*8)  with  a,  b,  and  0 functions 
of  r and  ph.  The  criterion  for  stability  becomes 

4(t-rt)~Z4r(i-r*)ph+(z7-48rt+  3fr*)phx+  4r($-5rl)pW+  /2r*p€h+)  0 (3 .2.3.10) 

Values  of  a,  b and  c are  listed  in  Table  3 .2,3=1  for  various  ph.  The  ratios  of 
Tg/T  and  Ag/A  are  also  plotted  against  ph  in  Figs.  19  and  20.  It  can  be  seen 
that  the  error  in  period  increases  with  time  interval  h and  the  damping  factor  r. 
The  amplitude  ratio  is  less  than  1 for  systems  with  slight  damping  and  greater 
t'-an  ] for  system  with  higher  damping  factor  r. 
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3.3  Numerical  Solution  of  Differential  Equations 

< 

?.3.1  Euler’s  and  Modified  Euler  Method 

Euler's  approximation  is  based  on  Vue  assumption  that  if  y is  expressed 
as  a function  of  x by  the  equation  ^ = f (*-)/)  , the  increment  in  y corre- 
sponding to  an  increment,  Ax,  in  x is  given  approximately  by  the  equation 
Ay  = /f*.  y)AX  , the  value  of  f(x,y)  being  that  at  the  beginning  of 
the  interval  Ax.  Applied  to  the  problem  of  free  vibration,  with  damping, 
governed  by  the  differential  equation  y + 2 rpy  t pxy  * 0 , the  formulas  for  new 


displacement  and  velocity  at  end  of  a time  interval  are  found  to  be 


y ■*  y*  + V* 


1.. 


(3. 3«i.l) 
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and  y,  = ( t-  2rph)  y6  - phy*  (3.3.1. 2) 

When  this  procedure  is  carried  on  in  a step-by-step  manner,  uair.g  the  displacement 
and  velocity  found  from  the  previous  time  interval  as  the  initial  condition  of  the 


new  mime  intervtu.,  tue  amei-Buto  equation  usCOiuc  = 

//v/  - 2(hrph)y*  t (t-  Zrph  + P'h')yn.,  **  « 

TL«  solution  of  this  equation  is  ^ 

y*  = G’(y>  + yj==3—Sm, yt) 

where  ^ - ^/^2rph  + p*hl 

and  /i  ® arc  <s/n  - arc  cos  - 


!-rph 


(?.3. 1-3) 


(3. 3. i.4) 


(3*3.1  5) 
(3* 3. 1.6) 


The  error  is  of  the  second  order  and  the  amplitude  damps  out  gradually 
with  time.  There  is  no  limiting  criterion  for  stability  in  case  of  free  vibration 
with  no  damping.  Any  time  interval  will  obtain  oscillatory  response.  For  free 
vibrations  with  damping,  the  criterion  for  critical  damping  is  the  same  in  the 
numerical  solution  as  in  exact  solution,  i.e,  aperiodic  motion  at  critical  damping 
occurs  when  r ■ 1. 

Figs.  21  and  22  show  the  errors  in  period  and  amplitude  at  pt  «*  1 for 
various  ph. 

The  modified,  method  of  Euler  takes  the  true  average  value  of  dy/dx  o'^er 
an  interval  instead  of  dy/dx  at  the  beginning  of  the  interval  for  the  equation 
& - fCx-y)  - ’Ibis  method  gives  a more  accurate  value  of  the  increment  of 
y due  to  increment  of  x than  Euler's  original  method  and  the  error  is  of  third 
order  of  (^xp.  This  method  may  be  represented  by  the  formula 

Af  K i(A'+  &*)  (3  3.1-7) 

where  ^ s 

- f(Xd  + Ax,  y0  + &y)Ay 


When  applied  to  the  problem  of  free  vibration  with  damping,  the  displacement  end 
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velocity  at  the  end  cf  an  Interval  now  become 

//*/  “ 0-  + (>-  rph)by„ 

yn+<  * [b 2rph - (i'irtjp'h1]'/*  * (-/  + rph)/?hyn 

This  p.vucBdure  may  be  represented  by  the  difference  equation 

ytH-tfhrfb-irrfftjy*  tl!'2rph+Zryh*-rpft+iph*]y„  * e 

the  aolution  of  which  is  found  to  be  . 

y*  « pn(ytu:n/4i 

where  R * 4/7-  2rpt?  + zr*y>'h*-  tfh'+iP'h* 

Again,  there  is  no  limiting  criterion  for  stability  in  the  case  of  free 
vibration  without  damping,  but  the  amplitude  of  vibration  seems  to  damp  out 
gradually  if  the  range  is  carried  too  far.  In  damped  vibration  problems,  the 
numerical  method  has  the  same  criterion  for  critical  damping  as  the  exact  solution. 

Pigs,  23  and  24  show  the  errors  in  period  and  amplitude  for  various  pb 
of  this  method. 

30»2  Runge,  Heun  and  Kutta's  Third  Order  Rule 

Various  formulas  have  been  devised  for  numerical  integration  of  the 
differential  equation  % ■ by  Runge,  Heun  and  Kutta.  These  oethode, 

accurate  to  the  third  order  of  Ax,  are  summarized  as  follows; 

1.  Runge' s Original  Formula? 

A/  a A*/+  (&*+  Aw)-  A""j  (5. 3. 2.1) 

a/  = f(*.y)-i ix 

A*  = f(*  + AX,  y+A')‘*x 

a"  = f(y+  ax,  y 

4"'=  f(>*  &*,  y*i*')-Ax 


(3* 3. 1.8) 
(3. 3- 1-9) 

(3.J.1.10) 

(3.3.I.U) 

(3-3.1.12) 

(3.3.1.13) 


where 


2. 


2.  Eeun; 

Ay  « $(*?*  3 a") 

(5. 3-2. 2) 

where 

a'  = f(x.y)-*x 

A*  a. 

f(x+fAX,y  + jA*)’A  X 

3.  Kutta'd  Third  Order  Rulet 

Ay  a (&'-  4A* + A**) 

(5- 3- 2. 3) 

where 

A'*  f(x.y)'&x 

A"*  f(y  + i&xt  y+’i[A')'AX 
f(*+AX,  y+2A*- Af)'A* 

All  these  formulas  when  applied  to  free  vibrations  with 

or  without  damp- 

log  yield  the  same  results.  When  applied  to  forced  vibrations , their  results  have 
alight  discrepanciss  but  all  contain  errors  of  the  fourth  order.  It  is  hard  to 
say  which  of  the  above  formulas  is  beBt  because  the  agreement  with  the  exact  solu- 
tion depends  on  the  type  of  forcing  function,  damping  coefficient  and  time  interval 
in  a very  complicated  manner. 

Considering  the  free  vibration  problem  with  damping,  the  following 
discussion  is  common  to  all  three  of  the  methods  above  mentioned. 

The  displacement  and  velocity  at  end  of  a time  interval  are  found  to  be 
/-0m  {p'h,ifp,¥)y  + { t-rph  - hye  (5. 5. 2. 4) 

These  equations,  if  applied  successively  by  using  the  final  value  of  previous 
step  as  the  initial  value  of  the  new  step,  lead  to  the  difference  equation 

>,/  - [t-Zrph-ti-lrtyhU  ■fC3-4f,JpW]y„ 


0. 3-2. 6) 


3$. 


The  solution  of  this  equation  is  , 

> ■»  &"( y. CtS n/  f <y ) (5-3.S.7) 

vfc,re  R * ./i^trph  1 2 ,‘p'h‘-  f H^b1- -jiO- &><■)?*-  £??+ fir*' 

?k v- rPh - to-^p^/nu  (5'5‘S-8) 

. tnu,  <-  - iO-‘r‘)t>V  + i(3-4rl)!>,h> 

* (3.3. 2.9) 


As  before,  these  methods  have  no  limiting  criteria  of  time  Interval  for 
stability  in  problems  of  free  vibration  without  damping,  nor  there  la  any  dis- 
crepancy in  the  damping  coefficient  r at  the  point  of  critical  damping. 


3.3.3  Kutta's  Fourth  Order  Rules 

The  formulas  vith  error  of  order  (Ax)^  derived  by  Kutta  are: 

1.  Kutta's  Simpson's  Rule 

Ay  " *"*)  (3. 3. 3.1) 

where  A*  • f /t)'AX 

&**/  (*•+ £**'  /•  * £*')'A* 

&*•  / (xt+i*xf  y*  + itf)'Ar 

A***  f (*,  + AX,  /0  f A")  'AX 

2.  Kutta's  Three -eighth Rule 

Ay  c f(*'+  Ja"  + 3a*'+  ***)  (3. 3- 3. 2) 

where  A*  * f AX 

a**  fCxt+j**,  y*  + i&O'A* 

£?•  f (Ve  + i' A* , y*  + AX 

A*9*  f(x6  + Ax,  y»+  A"'-  &"+  A'  )a  X 


lie  difference  in  results  between  these  two  rules  occurs  when  applied  to 
free  vibration  problems.  For  forced  vibrations  the  results  will  differ  slightly 
out  both  are  of  the  same  order  of  error.  The  agreement  with  the  exact  solution 
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varies  with  r,  ph  and  t. 


The  displacement  and  velocity  at  end  of  a time  interval  in  the  problem 

of  free  vibration  with  damping  are  found  by  either  of  the  two  rules  to  be 
V.  =r  [!-fpfy*{rpW+£(l-4r’Mi*ly. 

* [i-rph- {0-4  r*)f>*hl*  £(/-  2r*)p'h9]  hy0  (3.30*3) 

y,  « [-!  4 rph  4 {(i'4-r*)pW-  { (/-  2r*)fth3 ] p*hy0 

+ \l-2rph £(‘-f2r\/C>rjp4h4]y6  (J 

This  nay  be  put  in  difference  equation  form  for  the  step-by-step  method,  as  follows 

/*/ '/?- - (bzr*)^  iC3-4r*)pty+  jfO-  * <*)p4h*]yt 

+ [\-2tph+  2ryh9-  §r4p4h4+  £ (t-tr^pth* 

-■jiO-irtpW-  firpW*  jfcfiW]*.,  - o (300.5) 


(3000) 


The  solution  Is  found  to  be 


/* 3 Rn( y us*/'  + v*) 

vhera  R - [i-zrph  4 zrtpb*-  fr*pih3t§r4p4^4 

+ £(f-4r*y/,*-  jiO- br'Oy*- fop'll1  tJftpW]* 


A * 4/cstrt 


= &SCC0S 


/-  rph  - 


t-t-ph  - j 0-zr*)pVi*4  £(3-  Ar^jp*^  (t-  Sr'+Srfpj;* 


(300-6) 


(30.3.7) 


(3OO-0) 


Lite  all  other  methods  of  this  kind,  the  fourth  order  rules  have  no 
limiting  criteria  for  stability.  Any  time  interval  can  be  used  according  to  the 
accuracy  desired.  When  the  system  reaches  its  critical  damping,  i.e,,  r o 1,  the 
numerical  solution  becomes  aperiodic.  A difficulty  associated  with  the  method  is 
that  the  amplitude  will,  gradually  damp  itself  out  oven  in  the  case  of  an  undamped 
system,  and  the  method  is  therefore  not  desirable  for  a long  period  of  tine . 
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IV.  DISCUSSION 


4.1  Accuracy 

It  is  usually  convenient  to  compare  the  accuracy  or  various  methods  by 
the  order  of  time  interval  h involved  in  the  error  produced.  The  final  displace- 
ment over  a time  range  computed  by  numerical  methods  can  he  expanded  into  a ' 
polynomial  of  ph  and  then  compared  term  by  term  vith  a standard  pover  nerieo 
derived  from  the  exact  solution.  The  error  occurring  in  each  step  of  the  step-by- 
step  method  may  thus  be  observed,  especially  as  to  the  order  of  ph  involved. 
Considering  firBt  the  case  of  free  vibration,  with  damping,  which  is  governed  by 
the  differential  equation 

y * ?rpy  + p*y  * o (2.1.3) 


its  solution 


y - t^(y.  cos  JhT'pt , pt ) 


(2.1.4) 


may  be  expressed  in  a power  series  for  a time  interval  h 

X/  * X [ 4rl)p*h*~ £ (/- z r*)?#*  •“•] 

+ ynfi[i'rph-fc(t-4rt)jthl+  fc  ft-/t  f V /£  r4)^4*  ■ • • (4 .1.1) 

The  following  is  a collection  of  equations  which  show  the  accuracy  of  . 
different  numerical  methods: 

Constant  acceleration  method: 

* [i-t?W]y»  + [f-  rph]  hy„  (3.1.1.15) 

Newmark's  p-method: 


* [hiPtil+zprplhi-i(4firl-fi)ph4’ 

f i(4pr*-Qr-  4pr)p^s-  j (tfirtf/L  iprt+fifh V • Jy, 
+ f / - rph  + ( 4 fir1- zrfi)p*h* 

+ (tr4^- zr*$  - try  + /$%)p<h4+ 


(4.1.2) 


39- 


Fcr  p s 1/4  (Timoshenko's  modified  method-) 

3 ft- tpW  fifth* ]y„  + [i-  rph  - £( ’-trfp'h'*  --jhy,  (4.1.3) 

For  p a 1/6  (Linear  acceleration  method. ) 

r,  ±A(i.  r fiLi  , J_/.  *rt)*4Li. 1 ,/ 

» [i- + Jp"'  T ]Z[‘~  + r *r  “ r J'/ff 

* [/-rph-  t(/-*r,Jf,h,f  *f(*-t'l)pthif--‘]hy„  (4.1.4) 

For  p q 1/12  (Salvadori's  method.) 

/«-  = [/-  tpiu  ijfh1*  £(/-4r>Jp4t‘-  ■£(,.  3,i)ffhsf  ]y, 

For  p m 0 (Levy's  method.) 

“ [lmif^l]yn  f [!m  rph]hy*  (4.1.6) 

Parabolic  acceleration  method; 

/m  « [ I-  tpfri-fpVf&O-tr’tyV*  yffifa-Dirip’t**  ■ -Jy, 

[l-rpb-t(t-At')?h'+  i(i-2r'yP-’fi(l-lirli3lri)ith\.'.]t,}ll  (ltil  ?) 

Houbolt'8  method. 

>+/  ® [/-  4 ji(l-4r%)p*ii* f j y? 

* [/-rph- f(i-4rt)p*hti-  f (i-2r*)pW+  ■■■•]hy,  (4.1.8) 

Euler' a method; 

* hy»  (3.3.1 .1) 

Modified  Euler  Method. 

>W  - ( i-  )y«  t ( / - rph)!jyn  (3. 3. 1.6) 

Funge,  Heun  and  Kutta’s  third  order  rule; 

>/  = //- tfM'+j-fh*)#,  + [/- rph - iO-trtyh1] hy*  (3. 3-2. 4) 

Kutta's  method: 

/r)4i s [t~  i ftp*  jfpW  f f/~ 4r yp*h4] yn 

f [/- rph  - ’fc(t-4t't)plh'l+  i(f-2ri)p'ihi J hyt  (5-3. 3-3) 

From  the  above  listings  of  equations  for  various  methods,  one  may  observe 


1.  The  constant  acceleration,  linear  acceleration  and  parabolic  accel- 

03  4 

eratlon  methods  have,  respectively,  an  accuracy  of  the  order  of  b h and  h in 
both  displacement  and  velocity.  The  accuracy  of  Lsvmark's  S-method  depends  on 


1 


40. 


IV 


the  choice  of  0 . 

2.  Levy's  and  Houbolt's  methods  have,  respectively,  an  accuracy  of  the 
order  of  h and  hr  in  both  displacement  and  velocity  responses. 

5.  Salvador! 's  method  is  of  order  h^  only  in  displacement  response  and 

2 

only  for  undamped  systems . If  r and  y0  are  not  zero,  thi.9  method  is  only  of  h 
accuracy  unless  a good  Interpretation  of  yQ  is  made. 

4.  Euler's  method  is  of  the  first  order ; modified  Euler  method, -second 
order;  Kunge  and  Heun,  third  order;  and  JCutta's  Simpson's  rule  or  three-eighth 
rule,  fourth  order  generally. 

5.  In  some  special  cases  the  accuracy  of  the  above  methods  may  be  pro- 
moted one  more  order.  These  vill  be  listed  in  Chapter  5 . 

In  the  case  of  forced  vibrations,  the  error  vhich  enters  the  particular 
solution  also  governs  the  accuracy  of  the  method  since  it  determines  the  constants 
for  the  initial  conditions.  As  before,  Nevmark'c  linear  acceleration  method  is  of 
third  order  accuracy  for  any  system  vith  or  -without  damping.  For  an  undamped 
system  starting  at  rest,  Salvadori's  method  haB  an  accuracy  of  fifth  order.  Euler^ 
modified  Euler,  Runge-Heun  and  Kutta's  methods  are  Btill  of  first,  second,  third 
and  fourth  order  respectively. 

The  comparison  of  errors  by  polynomials  is  only  good  for  ph  less  than  1, 
leoetude  the  error  vould  othervise  be  dominated  by  the  higher  povers  of  ph  vhich 
would  vitiate  the  analysis.  Nevertheless,  usual  practice  indicates  that  a large 
time  interval,  say  from  l/4  to  l/2  of  the  natural  period  of  vibration  (ph=1.5708 
to  J.l4l6)  is  highly  desirable  in  rapid  and  rough  estimations.  The  effect  of 
large  time  intervals  is  shown  in  the  graphs  of  Appendix  J • 

4.2  Propagation  of  Errors 

The  preceding  article  concerns  the  accuracy  of  displacement  response  of 
various  methods  in  ono  single  step  of  operation,  tne  error  indicated  in  the 
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expanded  polynomial  being  compared  with  a common  initial  condition.  But  if  step- 
by-step  evaluations  are  set  up,  with  a common  initial  condition  for  the  firBt  step, 
the  initial  condition  of  the  second  step  will  contain  error  which  Ls  different  in 
various  method?  Aft*.'?  a chain  of  step?  completed,  the  error  propagates  in 
different  ways,  divergent  or  convergent,  accumulating  or  self -eliminating,  aa  & 
function  of  time,  length  of  time  Interval,  natural  frequency  and  damping  factor. 

Shown  in  previous  chapters  are  the  equations  for  yn  and  in  Appendix  3 
are  graphs  of  Tg/T  and  Ad/A  for  various  methods . These  may  serve  for  an  estimate 
of  the  propagation  of  errors 

Errors  are  of  two  fundamental  types,  error  in  period  and  error  in  ampli- 
tude. The  error  in  period  or  in  frequency  is  solely  due  to  the  discrepancy  in 
phase  angle  which  is  directly  proportional  to  the  lapse  of  time.  The  error  in 
amplitude,  on  the  other  hand,  chiefly  depending  on  the  factor  Rn  in  the  equations 
of  yQ,  is  an  exponential  function  of  time.  From  the  equations  for  yn  derived  and 
listed  In  Chapter  3,  it  is  evident  that  the  condition  for  no  error  in  period  ls 
ph  •»  n and  that  for  no  divergent  error  in  amplitude  is  Rerph  * 1.  All  pseudo 


periods  and  pseudo  amplitudes  have  been  compared  with  true  periods  and  true  ampli- 
tudes in  the  preceding  chapter.  The  ratio  of  periods  is  given  by  the  relation 

■ f (4.2.1) 

and  the  amplitudes  by  A0/A.  There  are  two  ways  of  comparing  the  amplitudes.  The 
first  one  concerns  the  ordinate  of  the  envelope  which  prescribes  the  periodic 
response,  while  the  second  way  takes  account  of  the  peak  amplitude  in  the  first 
cycle  of  vibration,  subject  to  a certain  initial  velocity. 

In  the  first  way  of  comparison,  one  finds  that 

*>  ■ ,,  0 


neglecting  the  magnifying  effect  cf  initial  velocity  in  some  cases. 


Since  the  ratio  of  amplitudes  is  ■'ot  constant  with  time,  it  is  reasonable 
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to  compare  the  methods  at  a certain  designated  time, 

Then 

Ml)  = oA,  tr 

Note  that  the  ratio.  A. /A  at  ct*l  here  dees 

o ■ 

with  Lhe  pseudo  amplitude  actually  computed  which  is 
factor  R,  hut  also  usually  by  the  change  ir.  velocity 
spense  depends  on  ph  as  stated  in  Art.  3 1.4,  and  sometimes  on  the  initial 

. t , 

dieplucement  as  in  £4.  (3*1. 3.18)  However,  when  the  time  interval  is  not  close 
to  the  criterion  of  stability,  this  error  in  amplitude  due  to  initial  velocity 
response  is  of  a constant  nature,  doing  much  less  harm  than  the  exponential 
factor  R after  a considerable  lapse  of  time,  and  is  not  taken  into  consideration. 
Therefore  the  ratio  Ag/A  for  the  envelopes  of  the  periodic  curves  is  still  useful 
for  judging  the  convergence  and  divergence  of  errors. 

Nevmark's  p-oethoa  for  all  values  of  p from  0 to  1/4,  together  with 
Timoshenko'3  modified,  Keva&rk's  linear  acceleration,  Levy'B  and  Salvador i's 
methods,  have  kjk  ■ 1 in  free  vibration  of  an  undamped  system.  Constant  and 
parabolic  acceleration  methods  are  the  ones  which  have  divergeng  amplitudes  while 
the  others,  including  Houbolt's,  Euler's,  modified  Euler,  Runge's,  Heun's  and 
Kutta's  methods,  have  convergent  aagplitude  although  the  rapidity  of  divergence  or 
convergence  is  different. 

The  presence  of  damping  may  add  complication  to  the  analysis.  Plottings 
of  Ag/A  at  pt  a 1 shown  in  Appendix  3 are  self-explanatory.  The  relation  between 
Ae/A  at  any  time  and  As/A  at  pt  “ 1 is  shown  iu  Fig.  29. 

In  the  second  method  of  comparison,  peak  amplitudes  in  the  first  cycle 
nf  vibration  due  to  an  initial  velocity  y0  a p are  compared  and  shown  in  Appendix 
3.  This  cop.tpariaon  may  be  of  more  interest  in  practical  problems  cf  vibration 
since  it  gives  actual  amplitudes  of  periodic  motion.  The  magnifying  effect  on 
the  sine  term  cf  the  general  equation  is  generally  taken  into  consideration  except 


say,  at  t a l/p,  i.e.  pt  0 1. 

(4.2.3) 

not  mean  the  comparison  made 
not  only  affected  by  the 
response.  The  velocity  re- 
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in  the  parabolic  acceleration  method  and  Eoubolt's  method  where  the  magnification 
factor  depends  also  on  the  treatment  of  initial  conditions.  Note  that  the  magni- 
tude of  the  peak  amplitude  does  not  depend  only  on  the  value  of  R,  but  also  on  the 
error  in  period  as  well. 

Error  in  period  is,  as  a rule,  constant  for  a given  ph.  Generally 
speaking,  the  ratio  T /T  increases  or  decreases  with  the  broadening  of  time  inter- 

O 

val.  Ts/T  } 1 indicates  a larger  pseudo  period  or  a retardation  of  phase  angle 
and  vice  versa. 

4 .3  Stability  and  Convergence 

The  applicability  of  different  available  techniques  places  some  limita- 
tions on  the  time  interval  used,  not  only  as  regards  accuracy,  but  also  as  regards 
stability  and  convergence.  All  acceleration  methodB  have  a limiting  criterion  of 
convergence  because  of  their  iterative  procedure.  All  acceleration  methods, 
except  Timoshenko's  modified  method,  also  have  a limiting  criterion  of  stability 
beyond  which  aperiodic  response  will  occur.  This  has  been  discussed  in  Art.  3.1,4. 
(See  Eq.  (5. 1.4. 7)  and  (3. 1.4. 8).)  Larger  values  of  (J  provide  a wider  range  of 
time  interval  for  stability,  but  a shorter  range  for  convergence.  On  the  other 
hand,  when  p ■ 0,  freedom  from  the  convergence  criterion  is  obtained  at  the  loss 
of  range  for  stability.  The  presence  of  damping  will  also  affect  both  criteria; 
the  greater  the  damping  factor  r,  the  shorter  the  range  of  time  interval  available. 

The  difference  equation  methods  have  generally  no  limiting  criterion 
for  convergence  because  of  the  nature  of  the  procedure . However,  a criterion  for 
stability  still  governs  those  techniques  which  have  been  discussed  in  previous 
chapters.  Levy's  and  Sulvadoi'i's  methods  have  the  same  criterion  for  stability 
an  that  of  Rewmark's  £5 -method  when  £ = 0 and  (3  = i/12  respectively-  The 
Houbolt  method  criterion  for  stability  has  been  given  in  Eq.  (3. 2. 3. 10). 

The  methods  of  numerical  integration  described  in  Art  3«3  have  the 
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advantage  of  avoiding  limiting  criteria  both  of  stability  and  convergence.  The 
choice  of  time  interval  nay  therefore  be  made  according  to  the  accuracy  desired. 
It  is  not  advisable  to  use  any  time  interval  greater  than  l/p,  that  is,  ph  > 1, 


utneiiTise  very  misleading  results  msy  be  obtained  because  of  th*  fact  that  higher 


powers  of  ph  will  dominate  the  solution. 


The  criterion  of  critical  damping  in  the  exact  solution  is  r * 1,  while 
only  some  of  the  numerical  solutions  have  the  same  criterion.  Those  which  have 


proper  criterion  of  critical  damping  are:  Timoshenko ' s modified  method  at  RevmariPe 

P-msthod  for  0 ■ 1/4,  Euler's  method,  modified  Euler  method,  Runge's,  Heun's  and 
Kutta's  methods.  In  moat  of  other  methods  critical  damping  occurs  even  when  r<^i. 
The  linear  acceleration  method,  Salvadori's  method,  Levy's  method  and  Revmark's 
0-method  (with  the  exception  of  0 * l/4)  sure  all  of  this  group.  The  criterion  of 
the  constant  acceleration  method  may  be  higher  or  lower  than  the  actual  criterion, 
depending  on  the  product  of  the  natural  frequency  of  the  structure  and  the  time 
interval  used.  Houbolt'a  method  generally  exhibits  periodic  response  for  all 
values  of  r and  ph  except  in  some  cases  when  r is  greater  than  0.94  and  ph  very 
email.  Fig.  30  illustrates  the  criteria  of  critical  damping  for  various  techni- 
ques. The  region  above  a curve  is  that  of  aperiodic  response,  while  the  region 
under  a curve  is  that  of  'periodic  response. 


4.4  Procedures  of  Operation 

Acceleration  methods  require  an  iterative  procedure  starting  from  an 
assumed  value  of  acceleration  and  arriving  at  a derived  acceleration  by  use  of 
the  equations  of  motion  until  o close  agreement  is  obtained  between  the  assumed 
and  derived  values . With  a proper  choice  of  time  interval,  three  of  four  trials 
for  each  interval  of  time  vill  usually  be  sufficient  to  reach  convergence  in  a 
mv.ltl-degree~of -freedom  system.  The  time  consumed  in  completing  a step  by  an 
electric  deck  computer  la  about  nine  minutes  for  a tvo-degree-of -freedom  system 


and  about  twenty  mlnutos.for  a five-degree-of -freedom  one  afttir  the  equations  are 
set  up  and  a routine  form  is  made  The  result  of  each  step  is  Belf -checking  except 
in  the  special  ease  of  g ■ 0 in  Nevmark's  (3-method.  Both  displacement  and  velocity 
rtay  Berve  an  a supplementary  help  for  giving  a clearer  picture  cf  motion  and  for 
checking 

Difference  equation  methods  arc  faster  because  displacements  are  directly 
obtained  from  the  difference  equations  and  no  extra  work  to  obtain  velocities  is 
needed.  However,  they  suffer  from  the  absence  of  self-checking  procedures  unless 
an  additional  device  is  provided.  For  problems  of  multi-degree-of -freedom  systems, 
Houbolt  suggested  a recurrence -matrix  solution  in  which  the  equations  of  motion 
are  expressed  in  a recurrence  matrix  equation  and  solved  by  inverting  the  matrix. 
Salvadofi  expressed  the  equations  of  motion  for  every  three  adjacent  masses  so 
that  each  equation  contains  only  three  unknowns  and  may  be  solved  by  relaxation, 
trial  and  error,  or  successive  approximations.  The  evaluation  of  displacements 
with  six  significant  figures  for  a five-degree-of -freedom  took  approximately  ten 
minutes  after  the  computations  had  been  standardized. ^ 

Hunge’s  and  Kutta  s methods  are  the  most  time-consuming  as  far  as  the 
use  of  an  electric  deek  computing  machine  is  concerned.  It  takes  more  than  thirty 
minutes  to  complete  a step  by  Kutta’ a fourth  order  formula  for  a tvo-degree-of- 
freedom  system.  Furthermore,  Binee  there  is  no  self-checking  of  calculations, 
mistakes  may  easily  be  introduced  into  the  computations  due  to  the  intricate  work 
of  cross-substitution  in  the  procedure . 
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The  general  results  of  this  study  are  tabulated  on  the  next  page  (Table 
5-1)  in  which  the  advantages  and  disadvantages  of  each  of  the  available  techniques 
are  listed.  Graphs  shoving  the  errors  in  period  and  amplitude  for  a range  of  time 
interval  from  0 to  about  half  of  the  natural  period  in  various  methods  are  given 
at  the  and  of  this  dissertation.  It  is  therefore  possible  to  choose  a suitable 
technique  for  a specific  problem  according  to  the  accuracy  and  amount  of  work  re- 
quired. In  general,  the  larger  the  time  Interval,  the  cruder  the  results  * Values 
of  ph  less  than  1 always  give  reliable  results  for  all  techniques,  but  variations 
will  be  great  when  ph^l,  and  these  graphs  may  be  found  useful  for  Judgement  vhen 
using  large  Intervale. 

In  ordinary  problems  of  vibratory  motion,  Nevmark's  p-method  is  most 
valuable  because  of  its  flexibility  in  application.  The  choice  of  time  interval 
may  be  made  for  the  desired  rate  of  convergence  and  accuracy  by  adjustment  of  the 
3 - parameter.  The  linear  acceleration  method,  a special  case  of  the  p-oethcd  for 

p»l/6,  la  most  consistent  in  degree  of  error  when  the  notion  1b  that  of  forced 

vibration  with  damping,  with  initial  displacement  and  velocity.  Timoshenko's  . 

method  is  test  applied  to  an  undamped  system  vhen  the  response  in  amplitude  is 

important.  The  constant  acceleration  method  and  Euler's  method  are  not  advisable 
owing  to  their  inaccuracy.  If  the  masses  in  motion  are  not  damped  and  have  no 
initial  velocity,  Salvador! 's  method  is  most  rapid  and  accurate.  JTor  rapid  and 
less  accurate  work,  Levy's  method  may  prove  useful,  but  care  should  be  taken  in 
the  treatment  of  initial  velocity.  Bunge  and  Kutta's  methods  are  noted  for  their 
accuracy  and  generality  in  application,  having  no  restrictions  with  respect  to 
stability  and  convergence,  but  they  are  handicapped  by  the  tedious  procedure  which 
is  not  generally  desirable  for  use  as  an  orainary  engineering  design  technique. 
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Table  5.1  --  Summarized  Result  of  Analysis 


Order  of  Accuracy 


Item 

No. 

Techniques 

Displace- 

ment 

Response 

Velocity 

Response 

Forced 

Vibra- 

tion 

Amplitude 
when  r»0 

<4.  • 

Constant  Acceleration 

2nd 

2nd 

divergent 

s. 

Tlao.  Mod.  Accel. 
Nevmark  p ■ l/4 

2nd 

(3rd  when 
r a o) 

2nd  2nd  constant 

(3rd  when  (Jrd  when 
r ■ 1/2)  r ■ 0) 

3. 

Linear  Acceleration 
Newmark  P * 1/6 

3rd 

(4th  when 
r - 1/2) 

3rd 

(4th  when 
r ■ 0) 

3rd 

constant 

• k; 

Parabolic  Acceleration 

4th 

4th 

4th 

divergent 

5' 

Nevmark ’s  p -method 

Depends  on  P 

constant 

6. 

Levy 

Revmark  P ■>  0 

2nd 

(3rd  when 
r - 0) 

2nd 

(3rd  when 
r o 1/2) 

2nd  constant 

(3rd  when 
r ■ 0) 

7. 

Salvador i 
Nevmark  p » 1/12 

2nd 

(5th  vii«n 
r ■ 0) 

2nd 

(3rd  when 
r - 1/2) 

2nd  constant 

(5th  when 
r ■ 0) 

8. 

Houbolt 

3rd 

(4th  when 
r - 1/2) 

3rd 

(4th  when 
r a 0) 

3rd 

conv. 

9. 

Euler 

1st 

1st 

(2nd  when 
r * 0) 

1st 

conv. 

10. 

Modified  Euler 

2nd 

(3rd  when 
r a 0) 

2nd 

2nd 

conv. 

U. 

Runge,  Heun 

3rd 

3rd 

(4th  when 
r a 0) 

3rd 

conv. 

12. 

Kutta  4th  Order 

4th 

4 th 

4th 

conv. 

(5th  when 
r - 0) 


Sumrarizad  Result  of  Analysis.  (Concluded) 


A 

i 


Table  9.1  - 
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Item 

No. 

:Ci iter ion  of 
Stability 

Criterion  or 
Convergence 

351?: — 
Check- 
ing 

fine' 

Conauap’ 

tien 

1. 

lq.  (3.1.1.13) 

no 

no 

less 

2. 

no 

rph*ifi)1  < 1 

yes 

fair 

3* 

Eq.  (J.l.3.13) 

tph+ifh't  / 

yes 

fair 

4. 

ph  < 3.357736 
for  r ■ 0 

ph<  2.020427 
fox  r 0 0 

yes 

more 

3- 

Eq.  (3.1.4.15) 

Eq.  (3.1.4.16) 

yeB 

fair 

6. 

ph<*(t-r*) 

no 

no 

less 

7. 

ph  < C(i- 

no 

no 

less 

0. 

Eq.  (3.2.3*10) 

no 

no 

Jair 

9* 

no 

no 

no 

less 

10. 

no 

no 

no 

fair 

11. 

no 

no 

no 

more 

12. 

no 

no 

no 

more 

f 


APPENDIX  1.  - - NOMENCLATURE 

The  letter  symbols  and  notations  used  in  this  thesis  are  defined  as 
follows  except  otherwise  noted, 

A o amplitude  of  the  moving  mass . 

As  ■ pseudo  amplitude  obtained  by  numerical  method, 

A,  B,  C,  » constants  determined  from  initial  conditions, 
a,  b,  c,  * coefficients  used  in  equations, 

e ■ 2.710....  base  of  'natural'  logarithms 

F ■ forcing  function, 

h « time  intervals 

k ■ constants . 

n « number  of  tine  intervals  taken  in  the  step-by-step  evaluation,  or 

subscript  to  designate  displacement  function  at  a particular  tine, 

88  y0,  yn  or  fQ,  *Q. 

p ■ circular  frequency  of  natural  vibration, 

r * factor  of  viscous  damping  in  terms  of  p. 

T ■ natural  period  of  a general  system,  in  general  T ■ Ji*/p. 

Tb  ■ pseudo  pe.iod  of  vibration  obtained  by  numerical  method, 

t ■ time . 

tn  ® time  corresponding  to  end  of  nth  time  interval, 

x - variable  in  general  equations . 

y •»  displacement  of  mass 

y m first  derivative  of  displacement  with  respect  to  time,  i.a.  velocity 

of  mass. 

y * second  derivative  of  displacement  with  respect  to  tics,  i.e,  accelera-. 

tion  of  mass. 

a exact  solution  of  differential  equation. 

* particular  solution  of  difference  equation  for  numerical  methods. 
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APPENDIX  1.  - - NOMENCLATURE  (Concluded) 

yp  a oxact  particular  solution  of  differential  equation. 

P ■ coefficient  measuring  proportion  of  acceleration  at  end  of  interval  in 

determination  of  displacement. 

1 1 » phase  angle  at  end  of  the  first  time  interval  obtained  by  numerical 

solution. 
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FIGURES 


Constant  acceleration  Method:  Displacement  ii 
Free  Vibration.  No  Damping 


onstaat  Acceleration  Kethod:  Error  in  Amplitude 


1.0 


Fig.  4 - Timoshenko's  Modified  Method  (or  Nevmark's  Method  for  p = l/4): 
Displacement  in  First  Cycle  due  to  y0  = p . Free  Vibration.  No  Damping. 


Timoshenko's  Modified  Method  (or  Nevmark’a  Method  for 

Error  in  Period. 


Modified  Method  (or  Newnark’s  Method  for 
Error  ia  Amplitude . 


Fig.  7 - Linear.  Acceleration  Method  (or  Nevmark's  Method  for  0 = l/6): 
Displaceitent  in  First  Cycle  due  to  yQ  = p.  Free  VibratioD.  No  IXamping. 


Stability  Criterion 


Linear  acceleration  Method  (or  Hewmark's  Method 
Error  in  Period. 


Fig.  10  - Parabolic  Acceleration  Method:  Displacement  in  the  First  Cycle  due  to 

yQ  =*  p-  Free  Vibration.  No  Damping. 


Effect  of  initic.l  conditions 
neglected.  Magnification  factor 
taken  as  unity.)  I | 


Parabolic  Acceleration  Method:  Error  in  amplitude 


A / A 


■npm 

■HUH 


Envelope  of  amplitude 
at  pt  ■ 1.  (Effect  of 
initial  condition  ne- 
glected. Magnificatia 
factor  taken  as  udit^ 
Peak  amplitude  in 
first  cycle  due  to 


Fig.  15  - Levy's  Method  or  Nevraark's  Method  for  0 o 0: 
Error  in  Amplitude . 
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Fig.  .16  - Nevmark's  g -Method  for  g * l/l2  and  Salvador! 's  Method: 

Error  in  Period 


Fig,  IT  - Nevmark's  p-Metbod  for  p <*  l/lS  and  Salvador's  Method: 

Error  in  Amplitude 


Houbolt's  Method  - Displacement  in  Firi 
Free  Vibration.  Ho  Damping 


A a/A 


mma 


It  m 0.73' 


r ■ 0 


Envelope  of  finiplitude  at 
pt  * 1 

Peak  amplitude  in  firBt  cycle 
due  to  y0  op. 


Fig.  26  - Range,  Heun  and  Kutta's  Third  Oi'der  Rule 
Error  in  amplitude 


If 


